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Abstract: This paper is aimed at separation treatment of low- and high-frequency components 
in polar motion forecasting and thenimproving time-series predictions. For the purpose, the 
empirical mode decomposition (EMD) is employed as a filter to extract low- and high-frequency 
signals from original pole coordinate data. The decomposition of the pole motion observations 
between 1986 and 2015 from the International Earth Rotation and Reference Systems Service 
(ERS) C04 seriesillustrates that the low-frequency fluctuations including inter-decadal, inter- 
annual, Chandler and annual wobbles and shorter-period high-frequency oscillationscan be 
separated from the observed time-series by the EMD. On the basis of separation, the least-squares 
(LS) extrapolation of models for annual and Chandler wobbles and for the linear trend are used for 
deterministic prediction of the low-frequency fluctuations, while the autoregressive (AR) 
technology is applied to forecasting the high-frequency oscillations plus LS fitting residuals. Pole 
coordinateforecasts are calculated as the sum of LS extrapolation and AR predictions 
(LS+AR).We have evaluated the accuracy of our long-term predictions (up to 1 year in the future) 
in comparison with the IERS official predictions in terms of year-by-year statistics of 5 years. It is 
shown that the accuracy of the LS+AR methodcan be significantly improved using a combination 
of the EMD and LS+AR (EMD+LS+AR). Also, the proposed prediction strategyoverall 
outperforms the IERS solutions. In addition, the predictions are compared with those from the 
Earth Orientation Parameters Prediction Comparison Campaign (EOP PCC). The comparison 
demonstrates that the developed schemeis a very accurate approach to predict polar motion. 
According to this study, it is concluded that polar motion predictions may be enhanced through 
separation treatment of different time-scale fluctuations and thus such processing seems to be 
necessary in pole coordinate prediction. 
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1. Introduction 
The Earth orientation parameters (EOP): universal time, , pole coordinates and 


nutation-precession corrections, describe the irregularities of the Earth’s rotation. 
Technically, they are the parameters which provide the rotation of the International 
Terrestrial Reference System (ITRS) to the International Celestial Reference System 
(ICRS) as a function of time. EOP derived from modern space geodetic technologies, 
e.g., Very Long Baseline Interferometry (VLBI), Satellite Laser Ranging (SLR) and 
Global Navigation Satellite System (GNSS), are not available for users in real-time 
applications resulting from the complexity of data processing. Since rapid and 
accurate EOP predictions are required for various fieldsassociated withreference 


frames such as interplanetary spacecraft tracking and navigation, positional 
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astronomy, geodesy and precise orbit determinations of artificial Earth satellites, it is 
essential to make high accuracy EOP predictions at least for a few days in the 
future''!.The currentpaper focuses on predicting, pole coordinatesup to 365 days in 
the future.Many different prediction approaches and techniques were used in the past 
to improve the prediction accuracy of polar motion. Most of these prediction 
approaches and techniques use only information within pole coordinate data, e.g., the 
combination of the least-squares (LS) extrapolationand autoregressive (AR) model 
(LS+AR) ? land combination of the LS extrapolationand neural network (LS+NN) 广 
5l. It has been reported that the LS+AR combination is one of the most accurate 
techniques for pole coordinateforecast [9. 

The basic problem in forecasting time-series is the necessity of separation 
treatment of high- and low-frequency variations'”!. This problem can be addressed by 
combination of the empirical mode decomposition (EMD)with the LS+AR model. 
Our work presents a hybrid EMD+LS+AR method for the enhancement of pole 
variation prediction. In this method, the low-frequency components consisting of the 
long-term trend, Chandler and annual wobbles determined by the EMD are first fitted 
and predicted by the LS extrapolation of thelinear polynomial and harmonic 
model,and then theshort-period and non-cyclichigh-frequency components together 
with the LS fitting residuals are modelled by the AR technique. The subsequently 
predicted high-frequency components and LS residuals areadded to the LS 
extrapolation in order to obtained the predicted values of pole coordinates.The 
comparison with the existing prediction methods and techniques shows that the 


presented EMD+LS+AR approach is atremendous tool for pole coordinate prediction. 


2. Prediction methodology 
2.1 Empirical mode decomposition 


The EMD can decompose a signal into a limited number of IMFs whose 
instantaneous amplitude and frequency are physically meaningful. It works on a 
simple assumption that, at any given time, the data may have many many coexisting 
simple oscillatory modes of remarkably different frequency, one superimposed on the 
other. Each component is defined as an IMF satisfying the following conditions™!. (1) 
In whole data set, the number of zero crossings and number of extrema must either be 
equal or differ at most by one. (2) At any data point, the mean value of the envelope 
defined using the localminima and the envelope defined using the localmaxima equals 


zero. With the above definition of an IMF, a signal can be decomposed into many 


IMFs through the following sifting process. 

Step 1: Identifyingall the local extrema in time-series . 

Step 2: Connecting all the local minima (maxima) by a cubic spline line to form 
the lower (upper) envelopes. 

Step 3: Calculating the mean value of the lower and upper envelopes, , and then 
computing the difference between and as a proto-IMF, , i.e., 

(1) 

Step 4: Examining whether satisfies the above-mentioned conditions. If it does, 
can be regarded as an IMF component. If not, the above sifting process is repeated. In 
the iterating process, is treated as the data in the next iteration: 

(2) 
After times of iterations, 
(3) 
and becomes the first IMF, i.e., 
(4) 
Step 5: Computing the residual and then treating itas the new data, i.e., 
(5) 

Step 6: Repeatedly applying the decomposition procedure to all subsequent . 

The procedure finally stops when the residual becomes a monotonic function or 
a function with only one extremum from which no more IMF can be extracted. 
Through the above procedure the original data are decomposed into IMFs and one 
residual obtained, , which can be either a constant or the adaptive trend, i.e., 

(6) 
2.2 Least-squares extrapolation 
In accordance with the EMD-based separated low-frequency components of pole 
variation data, a LS model may be fitted. As the sampling interval of the International 
Earth Rotation and Reference Systems Service (IERS)EOP C04 series are one solar 
day, we employ the following polynomial-sinusoidal function for LS fitting and 
extrapolation, i.e., 
(7) 
where the biasand driftof the linear trend, amplitudes, and phases , of the annual and 
Chandler wobbles have to be estimated by the LS method using the separated low- 


frequency fluctuation data, and, . 
3 


The model is applied to deriving the LS fitting residuals as for all used to fit 
the LS model. In the process of prediction, is extrapolated to obtain a deterministic 
prediction of the low-frequency fluctuations. 

2.3 Autoregressive model 
Stationary time-series may be modelled by the AR model of order [ARQ]. The 
technique allows one to capture the temporal dependence within stationary time- 
series. A stationary, zero-mean stochastic process (in our case ) is said to be an AR() 
if it satisfies 

(8) 
where is white noise with zero mean and variance , and are the AR coefficients. 

The order can be selected by the information criterion (AIC)"!. The AIC for an 

ARQ is based on the following statistics: 

(9) 
where is a likelihood function, is a vector of AR coefficients for an AR(), is an 
estimate of a variance and is a number of points of data. One selects which 
minimizes the AIC statistics. There are several approaches to estimate the AR 
coefficients; however, we use the Yule-Walker procedure. 

The fitted AR model to is applied to forecasting the LS fitting residuals plus the 
high-frequency oscillations of polar motion using the linear prediction operator. The 
forecasts for several steps in the future are computed recursively. 

2.4 Prediction strategy 

For prediction of , pole coordinates, an integrated prediction strategy referred to as 
EMD+LS+AR is applied as described in the flowchart of Figure. 1. The presented 
strategy comprises the following stages. In the first stage, the low-frequency 
components including the annual wobble, Chandler wobble and long-term trend in 
pole coordinates and the high-frequency signals containing the shorter-period and 
non-cyclic fluctuations extracted respectively from original data using the EMD 
algorithm as stated in Section 2.2. In the second stage, the extracted low-frequency 
oscillations including the annual and Chandler wobbles are fitted and then 
extrapolated by the above LS model. In the next step, the AR technique are used to 
model the LS fitting residuals plus the separated high-frequency components of pole 
variations, i.e., . The subsequently forecastedvalues of are then added to the LS 


extrapolation in order to obtain the predictions of pole coordinates. 


> LSmodel | +/ LS extrapolation 
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Figure. 1 Flowchart of the EMD+LS+ARstrategy for pole motion prediction 


3. Numerical examples and results 
3.1 Data description 


The IERS publishes daily values of EOP series with a latent of several days after the 
processing of observations of all advanced space geodetic techniques including 
GNSS, VLBI and SLR gathered at the permanent tracking stations distributed around 
the world. The EOP are estimated together with the station coordinates and velocities 
from each technique. The subsequently estimated EOP from all the above-mentioned 
techniques are then combined to derive the EOP C04 solutions". These finally 
combined EOP solutions are regularly published at the official IERS website. In this 
work, the , pole coordinate data from the IERS EOP 08 C04 series are used as data 
basis to carry out the numerical experiments. The present accuracy of , pole 
coordinates is about 30uas for the 08 C04 solutions. 

3.2 Separation of the low- and high-frequency oscillations in pole coordinates 
Taking the pole coordinate as an example, Figure. 2 shows the decomposed results of 
the 30-year pole variations from January 1, 1986 to December 31, 2015. It can be 
seen in Figure. 2 that the pole coordinate can be decomposed into a limited number of 
IMFs using the EMD algorithm, indicating that the oscillations of different time-scale 
in the component of polar motion are characterized by a small number of IMFs. As 
can be found in Figure.2, two of the modes found, namely and , show a noticeable 
seasonal oscillation in the range of a few hundreds of milliarcseconds, and thus can be 
considered as the two dominating oscillations along the and directions, namely 
annual wobble and Chandler wobble. In addition, we have found that the IMFs , and 
with longer periods and the residuals arising from the EMD analysis can be identified 
as the inter-annual and inter-decadal signals and the long-term “natural” trend in the 
observations, and the modes , and with the smaller amplitudes represent the shorter- 
period oscillations and the very high-frequency non-wobbling components. It is 
noteworthy that due to the mode mixing problem in the EMD the mode separation 


may be difficult when the signal contains components having adjacent periods, as it 
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occurs for the Chandler and annual wobbles. Since here we only focus on separation 
of the low- and high-frequency oscillations rather than extraction of a particular 
component, however, this mode mixing problem dose not limit our analysis. 

A spectral analysis of the IMFs plus based on fast Fourier transform (FFT) is 
exhibited in Figure. 3. In comparison, we also perform a similar analysis of the 
original data as given in Fig. 3. Visually, plus contains indeed contains the two 
dominant oscillating modes, i.e., Chandler and annual wobbles, with the same periods 
and amplitudes as the original observations do. Therefore, these finds confirm that the 
oscillating feature hidden in polar motion can be visually characterized by some 
oscillatory IMFs having particular physical meaning. To extract the high-frequency 
components of pole coordinates the IMFs , and are reconstructed, while the other 
modes are counted up to gain the low-frequency oscillations including the long-term 
trend, Chandler and annual wobbles.Figure.4 depicts the separated low-frequency 
time-series, high-frequency time-series of the , pole coordinates as well as the raw 
data. Their plots in Figure. 4 suggest that the low-frequency components including the 
Chandler wobble, annual wobble and long-term trend can be described by a low-order 
polynomial and harmonic model, while the AR technique is suitable for modelling the 
high-frequency terms. 
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Figure.2Original pole coordinate and its IMFs arising from the EMD algorithm 
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Figure. 3Spectral analysis of the oscillatory IMFs andoriginal pole coordinates 


(a) Original x, data (b) Original E data 


0.4 0.6 5 
© 02 ao 
o 四 
3 3 0.4} 
总 v 
三 0.0 3 
a a 
0.2 | 
E02 Š 
-0.4 1 1 i 1 1 0.0 1 1 1 L 1 
1985 1990 1995 2000 2005 2010 2015 2020 1985 1990 1995 2000 2005 2010 2015 2020 
(c) Low-frequency components of x, (d) Low-frequency components of H 
0.4 v T Tr + 0.6 r r r r 
8 02 g 
o S 
3 aoi) 
总 总 
3 00 2 
a a 
0.2 | 
E02 £ 
-0.4 H! 1 fi 1 | 0.0 . 1 1 n | | 
1985 1990 1995 2000 2005 2010 2015 2020 1985 1990 1995 2000 2005 2010 2015 2020 
(e) High-frequency components of x (f) High-frequency components of y 
|= P 
0.010 + 0.005 + r r r 
8 0.005 8 
3 3 
v v 
2 0.000 2 
a a 
£ -0.005 £ 
-0.010 n n 1 n 1 -0.005 | 1 n 1 L 1 L 
1985 1990 1995 2000 2005 2010 2015 2020 1985 1990 1995 2000 2005 2010 2015 2020 
year year 


Figure. 4 Original , pole coordinates, the low- and high-frequency components 
3.3Prediction results 

EOP 08 C04 series serve as data basis and the length of base data is 10 years. The 
predictions of , pole coordinates are generated by means of the two prediction 
schemes: LS+AR and EMD+LS+AR. The length of all predictions is equal to 365 
days.In order to compare with the official predictions from the IERS Bulletin A, 
whichcontains EOP predictions for 1 year into the future at daily intervals and is 


issued by the IERS Rapid Service/Prediction Centre at the U.S. Naval Observatory 


(USNO), all 365-day predictions of pole coordinate data are calculated every week at 
different starting prediction epochs from January 7, 2011 to December 31, 
2015because theIERS Bulletin A gives an advanced solution updated weekly.The 
comparison of the polar motion predictions is carried out by means of the analysisof: 
(1) the absolute values of differences between the pole coordinate data from the EOP 
08 C04 series and their forecasts, and (2) mean absolute error (MAE) prediction errors 
defined by the following! 
(10) 

where is the starting prediction epoch, isthe number of -step predictions, is the -step 
prediction of polar motion and is the corresponding observation of the IERS EOP 
C04 series. In our analysis, . 

The maximum absolute valuesof the differences between the pole coordinate 
data and their LS+AR predictions calculated at different starting 
epochsexceed98mas.The larger differences are for the selected polar motion forecasts 
with predictions lengths greater than 5 months (Figure. 5). In the case of the polar 
motion predictions published by the IERS Bulletin A, the maximum values of the 
differences between the data and IERS predictionsdo not exceed 91mas and those 
differences are smaller than those computed for our LS+AR scheme (Figure. 6). For 
the tested data span, the combination of the EMD and LS+AR model produces polar 
motion predictions with maximum absolute differences of less than 87masand the 
extreme perturbations are significantly smaller than for the LS+AR and IERS 
solutions (Figure. 7). 

The representative statistic MAE as well asthe corresponding error bars for our 
predictions in comparison with those forthe IERS official predictions are shown year- 
by-yearin Figure. 8. The results indicate that the accuracy of the LS+AR combination 
can be remarkably improved by incorporating the EMD into the prediction model, 
especially in long-term predictions.In contrast tothe IERSpredictions, our results 
obtained by the EMD+LS+AR approachappear to be comparable with those from the 
IERS for near-term prediction. However, for long-term prediction the proposed 
method noticeably outperforms the IERS solutions except individual prediction days, 
as can be seen apparently in Figure. 8. 

In order to furthermore demonstrate the effectiveness of the developed 


EMD+LS+AR prediction strategy, the results of thisstrategy arealso compared with 
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those from the EOP PCC lasting from October 1, 2005 to February 28, 2008, in which 
the prediction period and validation scheme had been clearly specified in advance. 
Under the same rules and conditions, the comparison with the contemporary 
prediction methods and techniques involved in the EOP PCC is shown in Figure. 9, in 
which the data time span is same, the statistics are referred to the EOP 05 C04 series 
and the MAE is used as the statistical measure. A list of participants and prediction 
methods can be found in [6]. 

From the comparison it can be said that the error of the predictions for 1 ~ 10 
days in the future by the EMD+LS+AR model is equal to theerror of the most 
accurate prediction method for forecasting polar motion, namely the LS extrapolation 
of the harmonic model and AR prediction (Figure. 9(a)) and Figure. 9(b)). For short- 
term (up to 30 days in the future) prediction, the prediction accuracy of the developed 
method is comparable with or slightly worse than the accuracy of the three most 
accurate prediction techniques used (Figure. 9(c)) and Figure. 9(d)). As to long-term 
forecast ,(up to 1 year in the future) the predictions are inferior to the most accurate 
results, in particular those of data. However, our results are significantly more 


accurate than the predictions contributed by the other participants (Figure. 9(e)) and 
Figure. 9(f)). 
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Figure. 5Absolute differences (in mas) between the pole coordinate data and their 


LS+AR predictions calculated at different starting epochs 
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Figure. 6Absolute differences (in mas) between the pole coordinate data and their 


TIERS predictions calculated at different starting epochs 
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Figure. 7Absolute differences (in mas) between the pole coordinate data and their 


EMD+LS+AR predictions calculated at different starting epochs 
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Figure. 8MAE for the LS+AR, EMD+LS+AR and IERS predictions of pole 
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coordinates from 1 to 365 days in the future as well as the corresponding error bars 
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Figure. 9MAE for the EMD+LS+AR predictions and those contributed to the EOP 


PCC for different prediction days 


11 


4. Summary and conclusion 
Our investigations show that pole coordinate time-series were decomposed into a 


small number of IMFs representing different time-scale fluctuations, e.g., seasonal, 
inter-annual and inter-decadal fluctuations. Therefore, the EMD can be considered as 
a filter to extract specific oscillatory signals in polar motion. For the purpose of 
separation treatment of different components in time-series prediction, the EMD is 
used as a filtering method for separating the low- and high-frequency components of 
pole variations in our work. The results have indicated that both the shorter-period 
high-frequency variations and the low-frequency oscillations including the long-term 
trend,Chandler and annual wobbles can be extracted through reconstruction of the 
partial IMFs. Based on the findings, the LS extrapolation of the low-order polynomial 
and harmonic model and AR technology are recommendedto predict the determined 
low-frequency and high-frequency components, respectively. The predictions of polar 
motion are generated by means of the combination of LS extrapolation and AR 
prediction. In order to evaluate the prediction accuracy of the proposed strategy, we 
have calculated the developed EMD+LS+AR as well as LS+AR predictions year-by- 
year at different starting prediction epochs betweenJanuary 7, 201land December 31, 
2015 in comparison with the IERS official forecasts. The figures of the MAE and 
absolute values of differences between pole coordinate data and their predictions 
illustrate that the presented hybrid approach indeedprovides the more accurate 
predictions more than the conventional LS+AR combination, implying that the 
performance of the LS+AR method can be enhanced by involving the EMD into 
prediction model.Our results gained by the proposed method are better than the 
predictions released by the IERS as well, although those for individual prediction 
days are exception. In addition, the comparison between our results and the forecasts 
from the EOP PCC confirms that the EMD+LS+AR combination is avery powerful 
toolto predict pole coordinate data. The main conclusionis that the prediction error can 
be decreased when the EMD is incorporated into the LS+AR combination. 
Furthermore, we conclude that the separation treatment of different time-scale 
oscillations may enhance polar motion predictions and thus such processing seems to 
be of necessity in time-series prediction. 
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利用 经 验 模 态 分 解 提高 极 移 预报 精度 
BAT, GM, BER 
(1. 中 国 科学 院 国 家 授时 中 心 ， 西 安 710600; 
2. 中 国 科 学 院 大 学 ， 北 京 100049; 
5. 中 国 科学 院 时 间 频 率 基准 重点 实验 室 , 西安 710600) 


摘要 : 经 验 模 态 分 解 (Empirical Mode Decomposition, EMD) 是 一 种 数据 驱 
动 的 自 适 应 非 线 性 、 非 平稳 信号 分 解 方 法 。 为 提高 极 移 预 报 精度 ， 将 EMD 应 
用 于 极 移 预 报 中 。 首 先 利 用 EMD 方法 对 极 移 序 列 进行 分 解 ， 获 得 极 移 的 高 频 
分 量 和 低频 分 量 ; 然后 采用 最 小 二 乘 (Least Squares, LS) 外 推 模 型 对 极 移 低 
频 分 量 进行 拟 合 ， 获 得 LS MERA; HARRAH HEI (Autoregressive, AR) 

模型 对 极 移 高 频 分 量 和 LS 拟 合 残 差 之 和 进行 建 模 预报 ;最 后 将 LS 模型 和 AR 
模型 外 推 值 相 加 获得 极 移 预 报 值 。 将 EMD 和 LS+AR 组 合 模 型 预报 结果 与 
LS+AR 模型 预报 以 及 地 球 地 球 定 同 参数 预报 比较 竞赛 (Earth Orientation 


Parameters Prediction Comparison Campaign，EOP PCC) 的 预报 结果 进行 比较 ， 
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结果 表明 ， 将 EMD 应 用 于 极 移 预 报 中 ， 可 以 明显 改善 极 移 预 报 精度 。 


关键 词 : 极 移 ; 预报 ; 经 验 模 态 分 解 ; 最 小 二 乘 ; 自 回归 模型 


